Criticality in Random Threshold Networks: Annealed Approximation and Beyond 
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Random Threshold Networks with sparse, asymmetric connections show complex dynamical behav- 
ior similar to Random Boolean Networks, with a transition from ordered to chaotic dynamics at a 
critical average connectivity K c . In this type of model - contrary to Boolean Networks - propagation 
of local perturbations (damage) depends on the in-degree of the sites. K c is determined analyti- 
cally, using an annealed approximation, and the results are confirmed by numerical simulations. It 
is shown that the statistical distributions of damage spreading near the percolation transition obey 
power-laws, and dynamical correlations between active network clusters become maximal. We inves- 
tigate the effect of local damage suppression at highly connected nodes for networks with scale-free 
in-degree distributions. Possible relations of our findings to properties of real-world networks, like 
robustness and non-trivial degree-distributions, are discussed. 



PACS numbers: 87.I8.Sn, 89.75.-k, 68.35.Rh 
I. INTRODUCTION 

Recently, research on complex dynamical networks has 
become increasingly popular in the statistical physicist's 
community The variety of fields where networks 

of many interacting units are found, e.g. gene regula- 
tion, neural systems, food webs, species relationships in 
biological evolution, economic interactions and the orga- 
nization of the internet, gives rise to the question of com- 
mon, underlying dynamical and structural principles, as 
well as to the study of simple model systems suitable as 
a theoretical framework for a wide class of complex sys- 
tems. Such a theoretical approach is provided, for exam- 
ple, by Random Boolean Networks (RBN) , originally in- 
troduced by Kauffman to model the dynamics of genetic 
regulatory networks in biological organisms Q,[5| . fn this 
class of models the state Ui € {0, 1} of a network node i 
is a logical function of the states of fcj other nodes cho- 
sen at random. The logical functions are also chosen at 
random, either with equal probability or with a suitably 
defined bias. A phase transition with respect to the aver- 
age number K of inputs per node is observed at a critical 
connectivity K c — 2 ^: for K < K c a "frozen" phase is 
found with short limit cycles and isolated islands of activ- 
ity, whereas for K > K c one finds a "chaotic" phase with 
limit cycles diverging exponentially with system size, and 
small perturbations ( "damage" ) propagating through the 
whole system. While a more detailed understanding of 
this transition by means of percolation theory [[| is still 
in its infancy, theoretical insight primarily was gained by 
application of the so-called annealed approximation in- 
troduced by Derrida and Pomeau M and successive ex- 



*This is the critical value for unbiased Boolean functions; a 
bias p in the choice of the Boolean functions, where p denotes 
the mean fraction of l's in the sites outputs, shifts the critical 
connectivity to K c — l/[2p(l — p)]. 



tensions ; the analytical study of damage spreading 

by means of this approximation allows for the exact de- 
termination of critical points in the limit of large system 
sizes N. 

Closely related to RBN are Random Threshold Net- 
works (RTN) , first studied as diluted, non-symmetric 
spin glasses [nil and diluted, asymmetric neural networks 
|^2| , ^3| . Networks of this kind also show complex dynam- 
ical behavior similar to RBN at a critical connectivity 
K c . For the case where the number of inputs per node K 
is constant, theoretical insight could be obtained in the 
framework of the annealed approximation |l4|]T5| ], but 
the extension of these results to more realistic szenarios, 
where the number of inputs per node is allowed to vary, 
fails because of the complexity of the analytical expres- 
sions involved. 

In this paper we undertake an alternative approach which 
allows us to calculate the critical value K c for RTN 
with discrete weights and zero threshold, but with non- 
constant number of inputs per node. Starting from a 
combinatorical investigation, we identify an additional 
degree of complexity in this type of RTN, which is 
not found in RBN: damage propagation depends on the 
in — degree of signal-receiving sites. The introduction 
of the average probability of damage propagation allows 
us to apply the annealed approximation and leads to the 
surprising result, that this class of random networks has 
a critical connectitivity K c < 2, contrary to RBN where 
one always has K c > 2. Numerical evidence is presented 
supporting our analytical results. Finally, we will give a 
short outlook on the complexity of phenomena in RTN 
near K Cl which are beyond the scope of the annealed 
approximation, and we will discuss how our findings 
could relate to properties of real-world networks, such 
as robustness ]l7|-p0|] and non-trivial degree-distributions 
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II. RANDOM THRESHOLD NETWORKS 

A Random Threshold Network (RTN) consists of N 
randomly interconnected binary sites (spins) with states 
it, = ±1. For each site i, its state at time t+1 is a function 
of the inputs it receives from other spins at time t: 

o- i (i + l)=sgu(/ i (t)) (1) 

with 

N 

fi(t)=J2 c w( t ) +h - w 

The N network sites are updated in parallel. In the fol- 
lowing discussion the threshold parameter h is set to zero. 
The interaction weights Cjj take discrete values cy = ±1, 
+1 or —1 with equal probability. If i does not receive 
signals from j, one has Cy = 0. 

The in-degree ki of site i thus is defined as the number 
of weights Cij with Cy ^ 0. If X denotes the average 
connectivity of the network, for large N the statistical 
distribution of in- and out-degrees follows a Poissonian: 

Prob(fc, = k) = e_i? ( 3 ) 
fc! 

This corresponds to the case where each weight has equal 
probability p — K/N to take a non-zero value. 

III. DAMAGE SPREADING IN RTN: 
DEPENDENCE ON THE IN-DEGREE K 

The most convenient way to distinguish the ordered 
(frozen) and the chaotic phase of a discrete dynamical 
network is to study the so-called damage spreading: for 
K < K c , a small local perturbation, e.g. changing the 
state of a single site cr, — > —Ci, vanishes, whereas above 
K c it percolates through the network. In RBN, the prob- 
ability p s that a site i propagates the damage when a sin- 
gle input j is changed (1— >0or0— >1) does not depend 
on the in-degree k of site i: if the possible 2 2 Boolean 
functions of the k inputs are chosen with equal probabil- 
ity, one has p s = 1/2; if the Boolean functions are chosen 
with a bias p, where p denotes the mean percentage of l's 
in the output of i, one has p s = 2p(l — p) @. Hence, in 
RBN damage spreading only depends on the out-degree 
of the perturbed site. In RTN, the situation turns out 
to be more complex: here we find that damage spread- 
ing strongly depends on the in-degree of the sites. In 
the following, we will use combinatorical considerations 
to derive the exact distribution p s (k) for RTN. 

Consider a site i having k arbitrary input spins, k £ 
{0, 1, 2, N}. Let k + denote the number of spins equal 
to +1, fc_ the number of spins equal to —1, hence k = 
k + + fc_ . The state cr, (t + 1) of site i at time t + 1 is given 



by eqn. (1) and (2). Let us first calculate the probability 
p s (k) that a change of the sign of one arbitrary input 
spin at time t reverses the sign of i's output at t + 1, i.e. 
that it leads to a-(t + 1) = -a. t (t + 1). 
The number of possible configurations of k spins is 2 fc ; 
as in each of these configurations k spins can be flipped 
(reversed in sign) , the total number of possible spin- flips 
is 

Ztotal =k-2 k . (4) 

Thus p s is defined as the number of spin reversals leading 
to a[{t + 1) = -<Ti(t + 1), devided by Z tota i- 
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FIG. 1. Combinatorical derivation of p a (k), here for k = 3. 
Reversal of the thick- lined inputs dj(7j(t) leads to a spin- flip 
of site i, i.e. o[{t + 1) — — cr»(t + 1) (thick- lined vectors). 
The total number of input configurations is 8 (see the right 
column showing the corresponding numbers of spin configu- 
rations), hence the number of possible spin-flips is 3 • 8 = 24. 
The total number of input-spin reversals leading to damage 
propagation is 2- (3 + 3) = 12. Thus we findp s (3) = 1/2. The 
generalization for k = 5, 7, 9, ... is straight-forward. 

For fe = 1 it is easy to see that p 8 = 1, whereas for 
k = 2 one gets p s = 1/2. For k > 3 we have to analyze 
even and odd k seperately: 

Odd k: By the definition of the transition function for 
site i (eqn. (2)) one recognizes that there are only two 
configurations in which a flip of a single input spin at 
time t can lead to a spin-flip of site i at time t + 1; these 
configurations are given by the condition: 

|A+-fc_| = l. (5) 

If k + = k- + 1 (i.e. <x,(t + 1) = 1) , the reversal of 
a positive input spin always leads to 0^(4 + 1) = — 1, 
whereas flips of negative input spins do not change the 
state of site i (for fc + = — 1 it is vice versa). Thus 
in both configurations (k + l)/2 spin flips of k possible 
lead to a[{t + 1) = — <7j(i + 1). The total number of 
configurations fulfilling eqn. 5 is 
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Hence the number of spin-flips leading to damage spread- 
ing at site i is Z &dfUp 1 (k + l)/2 and we obtain for 
fc — 3 , 5 , 7, ... 



n ... _ z oddflip ■ (fc + 1) _ {k + 1} V(fc + lyy 
Ps{k} - 2-Z total k~2 k (7) 

Fig. [l] demonstrates the above derivation at the example 
fc = 3. p s {k) for odd k is shown in Fig. ||. 
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FIG. 2. Local probability for damage propagation, p s {k), 
as a function of the in-degree k of the signal-receiving site, 
for odd k. The squares correspond to the exact result of eqn. 
(7), the dashed curve shows the Stirling approximation of eqn. 
(13). 

Even k: For even k, there are also only two configu- 
rations in which changing a single input of site i at time 
t leads to a[(t + 1) = —<Ji{t + 1), given by the conditions 



k + = fc_ (*) or fc_ = k + + 2 



(8) 



In the case fc+ = fc_ + 2 neither the reversal of a posi- 
tive nor the reversal of a negative input-spin changes the 
output of i, due to the signum function in eqn. (2). In 
the case (*), k/2 spin- flips (+1 — » —1) change the output 
of site i (+1 — * —1), in the case (**) (k + 2)/2 spin-flips 
(— 1 — * +1) do so. The number of configurations fulfilling 
eqn. (8) is given by 



Zevenflip = Z{*) + Z(**) 



/fc 

So we finally obtain for fc — 2,4, 6, .. 

k-Z(*) + (k + 2)-Z(**) 



k 

(fc + 2)/2 



Ps(k) 



2 • Z 



total 



k ■ 



k/2 



+ (fc + 2) 



(fc + 2)/2 



k ■ 2 k + 1 



= P.(k + l), 



(9) 

(10) 

(11) 
(12) 



as can be seen by some simple algebra. Using the Stir- 
ling formula for fc!, one recognizes that the asymptotic 
behavior of p s (k) for odd k is given by 



Ps{k) 



1 



k + 1 

Thus it turns out that 



lim lim p s (k) = 0. 



!'7T < 13 » 



(14) 



However, this does not mean that there is no damage 
spreading for k —> N: in each timestep, damage increases 
oc k ■ p s {k) ~ \fk ~ v^V" in this limit, which ensures the 
existence of a "chaotic" regime. 



IV. AVERAGE PROBABILITY FOR DAMAGE 
SPREADING 

Choosing an arbitrary network site i with at least one 
input and changing the sign of one input-spin will lead 
to a spin- flip of site i with probability p s (k), dependent 
on the number of inputs k. Now we are interested in 
the expectation value of this probability, i.e. the average 
probability for damge spreading (p s )(K) for a given av- 
erage network connectivity K. For large N, this problem 
is equivalent to connecting a new site j with state aj to 
an arbitrary site j of a network of size N — 1 with aver- 
age connectivity K; thus ki is increased by one. Hence, 
changing the sign of Oj at time t will lead to a different 
output of site i at time t + 1 with probability p s {ki + 1). 
Using the Poisson approximation and taking the thermo- 
dynamic limit N — > ao, this leads to 



(Ps)(K) 



-K 



fe=0 



fc! 



(15) 



Splitting the sum for even and odd k yields 



+ g(^- p * (2 ' +1) I 



Using the relation p s (k) = p s (k+l) for even fc one obtains 



(p.) (i?) =e^-i+^ ^ 2i -v s (2i + 1) 

l i=l 

1 i? 



,(2* -1)1 (2i)! 
Inserting (7) into this equation finally yields 



(17) 
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(Ps)(K) = e 



■ R l+-r — (- 



K 



(18) 



Formula (18) is the central result of this paper, as it 
allows for an analytic calculation of the critical connec- 
tivity K c of RTN. For practical use it is worth to notice 
that the sum in eqn. (18) converges very fast, setting 
the upper limit to ten is sufficient for a relative error 
A(p s )/(p s ) < 0(cor A ). 
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FIG. 3. Time evolution of the relative Hamming distance 
yt of two spin configurations with yo = 0.3 for different net- 
work sizes N with K = 1.7 (numerical simulations, ensemble 
statistics over 10000 networks / 400 networks (N = 32768)). 
For large N, one finds convergence against the prediction of 
the annealed approximation (solid curve). The inset shows 
the same for overcritical networks (K = 2.0). 



V. CALCULATION OF K c BY DERRIDA'S 
ANNEALED APPROXIMATION 

The still most powerful analytical approach to describe 
damage spreading in random networks of automata - 
and thus to calculate the critical value K c of K - is the 
so-called annealed approximation introduced by Derrida 
and Pomeau [Q. This approximation, originally applied 
to Kauffman networks (Boolean networks with constant 
number K of inputs per node), neglects the fact that 
the (Boolean) functions and the interactions between the 
spins are quenched, i.e. constant over time, and instead 
randomly reassigns inputs and functions to all spins at 
each time step. As a central result of this approximation 
Derrida and Pomeau derived a recursion formula which 
describes the time evolution of the Hamming distance 
<i((7, a 1 ) of two spin configurations a and a' . The normal- 



ized distance y = d(a,a')/N at time t + 1 for p s = 1/2 
and in-degree k is given by the recursion 



Vt+i 



(19) 



They also give a straight-forward generalization to net- 
works where k is not constant: 



Vt+i 



k 



1- (1-2/0* 



(20) 



where pk is the probability that a network site has k 
inputs. In the case of an average probability for damage 
propagation {p s ){K) depending on K, as for the RTN 
discussed in this article, this generalizes to 



+i ={ Ps )(K) Pk •[!-(! - 



yt) 



(21) 



In the following, we focus our discussion on random net- 
works, where each link has equal probability to take a 
non-zero value (as first introduced by Erdos and Renyi 
Jl6||). Thus, for average connectivities K <C N, pk is 
given by a Poissonian, leading to 



Vt+i = (Ps)(K)e 



-K 



yt) 



[K(l-yt)) h 
M 



(22) 

Fig. compares the time evolution of the Hamming dis- 
tance, measured in numerical simulations of RTN of dif- 
ferent size N, to the prediction of eqn. (22). 
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FIG. 4. The normalized distance y for t — > oo as a function 
of K, calculated by numerical solution of eqn. (23). One finds 
y(K) > for K > K c = 1.849. 
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FIG. 5. Average Hamming distance d t for two configura- 
tions a and <?' differing in one bit at time t — 1 as a function 
of the average connectivity (ensemble average over 10000 ran- 
dom networks with N — 128 for each data point). One finds 
excellent agreement with the prediction of the annealed ap- 
proximation (solid curve, d t , a „„ e aied = (p s )(K) ■ K). 

As well for subcritical as for supercritical K one finds 
convergence against the annealed approximation, but for 
K > K c the convergence is quite slow (see section 6). 
In the limit t — > oo we expect that the normalized dis- 
tance y(K) evolves to a constant value, i.e. yt+i — yt — 
y. Inserting this condition in eqn. (22) leads to the fixed 
point equation 



f(y) = y- ( Ps )(K) (1 - cxp [-K ■ y}) = 0. 



(23) 



Solutions y* of eqn. (23) for some values of the average 
connectivity K are shown in Fig. 4; evidently, there ex- 
ists a critical connectivity K c above which the Hamming 
distance of two initially nearby trajectories increases to 
a non-zero value, indicating damage spreading through 
the network. The exact value K c can be obtained easily 
by a linear stability analysis of eqn. (23). j/q = is al- 
ways a solution of (23), but is attractive (stable) only for 
K <K C . 

This fixed point is stable only if 



lim 



df(y) 



dy 



(y* + £)<!• 



(24) 



One has 



df(y) 



dy 



(y* +e) = (p.) (K)-K- exp [-K ■ s], (25) 



i.e. 



lim 

e^0 



df(y) 



dy 



(y* + e) = (p s )(K).K. 



Inserting this result into (24) we find that y$ 
attractive if 

{p s )(K)-K< 1. 



(26) 



is 



(27) 



Thus the critical connectivity K c can be obtained by solv- 
ing the equation 



(p s )(K c )-K c = l, 



(28) 



{p s ){K) given by eqn. (18). 

Eqn. (27) and (28) have a simple interpretation in terms 
of damage spreading: (p s )(K)-K is the expectation value 
(dt) of the Hamming distance - i.e. of the damage - at 
time t, if at time t — 1 a single spin is reversed, as this 
spin on average has K outputs to other spins, which all 
will propagate damage with probability (p s )(K). Fig. || 
shows that the prediction (dt annealed) of our annealed 
approximation agrees perfectly well with with the aver- 
age dt measured between two spin configurations in RTN 
with dt-\ = 1 (ensemble statistics). 
A numerical solution of eqn. (28) for the RTN discussed 
here yields the critical value 



K r = 1.849 ±0.001. 



(29) 



Thus we get the remarkable result that in RTN with 
h = and discrete interaction weights, marginal damage 
spreading - i.e. the percolation transition from frozen to 
chaotic dynamics - is found below K = 2, in contrast to 
RBN where one always has K c > 2. 



VI. BEYOND THE ANNEALED 
APPROXIMATION: COMPLEXITY IN RTN 

So far we presented numerical evidence that the 
annealed approximation correctly predicts the average 
damage spreading behavior in RTN and thus the crit- 
ical connectivity K c in the limit of large system sizes, 
however, this coarse-grained approach of course does not 
capture the whole complexity of the network dynam- 
ics near criticality. Concerning the statistics of dam- 
age spreading, even for quite large N scale-free distri- 
butions are found in a certain range around K c ; skewed, 
super-critical distributions, with an increasing maximum 
moving towards N/2 with increasing K, are found for 
K > 2.1 (JV = 8192, Fig. |). In the ordered regime, 
the distributions decay exponentially. Presumably, the 
good convergence of the average Hamming distance mea- 
sured in numerical simulations for K < K c against the 
annealed approximation directly reflects averaging over 
these exponential distributions with well-defined charac- 
teristic scale, whereas averaging over the the scale-free 
damage-distributions around K c leads to the observed 
weak finite-size scaling oc 1/ log (N) against the analyt- 
ical result (Fig. ||) at the order-chaos transition. Fur- 
thermore, for finite JV, even the transition from critical 
(scale-free) to supercritical distributions (Gaussian dis- 
tributions in the limit K — > JV) is rather smooth: even 
deep in the "chaotic" regime (e.g. K = 2.5, Fig. |^) 
damage becomes zero for more than 60% of the initial 
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conditions/ the networks tested, and the damage distri- 
butions are clearly bimodal. 

Attractor periods of RTN are found to be power-law 
distributed in a certain range of connectivities around 
K c , due to dynamical correlations between network sites, 
which are neglected completely in the framework of Der- 
rida's annealed approximation. We shall briefly discuss 
this: The average correlation Corr(i,j) of a pair of 
sites is defined as the average over the products ai(t)aj (t) 
in dynamical network evolution between two distinct 
points of time T\ and T^'- 



1 2 

Corr(i,j) = Yl ^ViW 



(30) 



T2-T1 

The global average correlation can be defined as 

^ N N 

( Corr ^ = N(N-1) S 5 |Corr(i, j) | with i ^ j, 



=1 j=i 



(31) 



neglecting trivial auto-correlations. 
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FIG. 6. Statistical distributions p(d) of the Hamming 
distances d(a, a) for two spin configurations a, a with 
d(t = 0) = 0.3 ■ 2V. Only the tails of the distributions are 
shown in this log-log-plot. Ensemble statistics was taken af- 
ter 100 dynamical updates over 5000 RTN with size N = 8192 
and Poissonian distributed in- and out-degree, testing 20 dif- 
ferent initial conditions for each network. In the ordered 
regime (K — 1) damage is suppressed exponentially. Near 
criticality and slightly above (K = 1.85, K = 2.0) the 
tails aproximately obey power-laws (log-binned data). For 
2.1 < K < 2.6, the distributions are bimodal (shown here for 
K — 2.5), with an increasing maximium at large damage val- 
ues. For larger K, the tails of the distributions become Gaus- 
sian, with a maximum approaching d — N/2 with increasing 
K. Not visible in this plot is the pronounced maximum of 
all distributions at d — 0: even for K — 3.0 damage finally 
becomes zero for about 56% of the networks/ the initial con- 
ditions tested; this maximum vanishes for K — *• N. 



In the ordered (frozen) regime, there are few, isolated 
(non-frozen) clusters with dynamical activity, but their 
dynamics is not correlated. A good measure here is 
the average correlation of non- correlated sites (Corr) nc , 
i.e. the average correlation of pairs of sites with < 
Corr(i, j) < 1 (Fig. [7]). (Corr) nc is almost zero in the 
ordered regime. Due to this uncorrelated dynamics of 
a few "active islands" in a frozen network the distribu- 
tions of attractor periods P decay oc exp (— P). Near 
K c , however, (Corr) nc shows a pronounced maximum, 
i.e. the dynamics of active clusters becomes strongly 
correlated; thus, at this percolation transition, attrac- 
tor periods show scale-free distributions. In the chaotic 
regime, damage spreading destroys most of the correla- 
tions, consequently, (Corr)„ c decays once again. 
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FIG. 7. Global average correlation (Corr} 3 (lines) and av- 
erage correlation of non-correlated pairs (Corr) nc (crosses or 
points) as a function of the average connectivity K for four 
different system sizes N. For each value of K sampled en- 
semble-averages were taken over 1000 RTN. For each network 
the average was taken over a subset of 5000 randomly chosen 
pairs of sites. Individual runs were limited to T ma x = 20000. 



VII. DISCUSSION 

We investigated damage spreading in Random Thresh- 
old Networks (RTN) with zero threshold and discrete 
weights. This kind of discrete dynamical network shows 
complex dynamics similar to Boolean networks, with a 
transition from ordered to chaotic dynamics at a criti- 
cal average connectivity K c . Using combinatorial meth- 
ods, the exact distribution p s (k) for local damage prop- 
agation was derived. It was shown that Derrida's an- 
nealed approximation can be applied to this class of mod- 
els; this theoretical analysis yielded the surprising result 
K c = 1.849, in contrast to RBN, where K c > 2. The 
ansatz proposed in this paper could offer a road-map 
for an analytical treatment of similar systems with addi- 
tional complexity (non-discrete weights, h ^ 0). 



G 



An interesting result of our studies is that in dynamical 
networks with sparse asymmetric interactions and some 
kind of threshold characteristics governing dynamics, a 
high number k of inputs per node can stabilize dynamics, 
as the effect of single, local errors is reduced like 1/yk; 
this, however, is counteracted by the overall topologi- 
cal randomness in the network, which allows for ordered 
dynamics only at small average connectivities K. This 
effect of "local damage suppression" at nodes with high 
in-degree ist demonstrated in Fig. 8, directly comparing 
damage distributions in RTN and RBN: damage spread- 
ing was investigated for networks with "flat" , scale-free 
in-degree distributions oc , wheras the out-degree fol- 
lows a Poissonian. For very flat in-degree distributions 
(a = 1), both RTN and RBN have a hierarchical, clus- 
tered structure and show ordered dynamics, but in RTN 
damage is stronger suppressed than in RBN, due to local 
damage suppression at highly connected nodes. 
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FIG. 8. Statistical distributions p(d) of the Hamming 
distances d(a, <x) for two spin configurations a, a with 
d(t = 0) = 0.3 • N. Only the tails of the distributions are 
shown in this log-log-plot. Ensemble statistics was taken af- 
ter 100 dynamical updates over 10000 RTN or RBN, respec- 
tively, with size N = 2048 and average connectivity K — 2.5, 
Poissonian distributed out-degree and in-degree-distribution 
p(ki n ) oc k~^" . The choice a = 1 leads to clustered, hierar- 
chical networks with ordered dynamics, but in RTN damage 
spreading is stronger suppressed: p(d) shows a faster expo- 
nential decay than in RBN. This is a direct effect of the local 
damage suppression due to highly connected nodes in RTN. 
For a = 2.2, the opposite is observed: Whereas RBN show an 
only slightly overcritical p(d) with a power-law regime prevail- 
ing over almost two decades, RTN show a strongly skewed, 
overcritical p(d). In this regime the local damage suppres- 
sion effect oc \Jkin is not strong enough to counteract chaotic 
dynamics increasing oc kout- 

For steeper in-degree distributions (a — 2.2 in Fig. 
8) this effect is not strong enough, and the dynamics 



in RTN is even "more chaotic" in RTN than in RBN. 
Concerning models including self- organization of net- 
work topology this could have important consequences. 
Topological evolution of RTN by a local coupling of 
control- and order parameters (connectivity and magne- 
tization/correlations of network sites) was introduced in 
p3j. Recent, more detailed studies of these models show 
that the self-organization processes, balancing the net- 
works in a regime of complex, non-chaotic dynamics near 
criticality, indeed favor highly-connected nodes to some 
extent, leading to deviations of the statistical distribution 
of in-degrees from a Poissonian H). Many networks in 
nature are characterized by non-Poissonian degree distri- 
butions (namely power-laws). Whereas for fast-growing 
networks like the internet models based on preferential 
linking provide explanations for the observed distribu- 
tions pjJ, for other dynamical networks - e.g. protein 
networks ]2lJ or gene networks ||| - convincing models 
are still not available; szenarios based on preferential link- 
ing have to make detailed a-priori-assumptions about the 
evolutionary processes leading to the observed structures, 
which in the case of biological networks usually cannot be 
falsified. Approaches based on network dynamics |^0|,^3| , 
setting the focus on the robustness of network-dynamics 
and -evolution, could provide more realistic szenarios, as 
some threshold characteristics usually is an intrinsic dy- 
namical feature of these networks. We expect that future 
research will underline the significance of RTN as simple 
"toy systems" yet able to capture the essentials of natu- 
ral dynamical networks: evolution of high robustness and 
non-trivial randomness. 
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